{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Example: Obtaining results from the eq and profiles objects\n", "\n", "Here we will take a look at how to access different results stored (or calculated using methods) in the `eq` and `profiles` objects. \n", "\n", "To do this we first need to run a static forward simulation in the MAST-U-like tokamak using previously saved coil currents to generate our results. Of course, the same methods below can be applied after an inverse solve!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Static forward simulation" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "import numpy as np\n", "\n", "# build machine\n", "from freegsnke import build_machine\n", "tokamak = build_machine.tokamak(\n", " active_coils_path=f\"../machine_configs/MAST-U/MAST-U_like_active_coils.pickle\",\n", " passive_coils_path=f\"../machine_configs/MAST-U/MAST-U_like_passive_coils.pickle\",\n", " limiter_path=f\"../machine_configs/MAST-U/MAST-U_like_limiter.pickle\",\n", " wall_path=f\"../machine_configs/MAST-U/MAST-U_like_wall.pickle\",\n", ")\n", "\n", "# initialise equilibrium object\n", "from freegsnke import equilibrium_update\n", "eq = equilibrium_update.Equilibrium(\n", " tokamak=tokamak,\n", " Rmin=0.1, Rmax=2.0, # radial range\n", " Zmin=-2.2, Zmax=2.2, # vertical range\n", " nx=65, # number of grid points in the radial direction (needs to be of the form (2**n + 1) with n being an integer)\n", " ny=129, # number of grid points in the vertical direction (needs to be of the form (2**n + 1) with n being an integer)\n", " # psi=plasma_psi\n", ") \n", "\n", "# initialise profile object\n", "from freegsnke.jtor_update import ConstrainPaxisIp\n", "profiles = ConstrainPaxisIp(\n", " eq=eq,\n", " paxis=8.1e3,\n", " Ip=6.2e5,\n", " fvac=0.5,\n", " alpha_m=1.8,\n", " alpha_n=1.2\n", ")\n", "\n", "# initialise solver\n", "from freegsnke import GSstaticsolver\n", "GSStaticSolver = GSstaticsolver.NKGSsolver(eq) \n", "\n", "# set coil currents\n", "import pickle\n", "with open('data/simple_diverted_currents_PaxisIp.pk', 'rb') as f:\n", " current_values = pickle.load(f)\n", "\n", "for key in current_values.keys():\n", " eq.tokamak[key].current = current_values[key]\n", "eq.tokamak[\"P6\"].current += 100\n", "\n", "# carry out forward solve\n", "GSStaticSolver.solve(eq=eq, \n", " profiles=profiles, \n", " constrain=None, \n", " target_relative_tolerance=1e-9)\n", "\n", "# plot the resulting equilbrium\n", "fig1, ax1 = plt.subplots(1, 1, figsize=(4, 8), dpi=80)\n", "ax1.grid(True, which='both')\n", "eq.plot(axis=ax1, show=False)\n", "eq.tokamak.plot(axis=ax1, show=False)\n", "ax1.set_xlim(0.1, 2.15)\n", "ax1.set_ylim(-2.25, 2.25)\n", "plt.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Dashboard of results\n", "\n", "It is worth playing around with the fields/methods in the `eq` and `profiles` objects yourself to see which quantities can be calculated from the plasma equilibrium using built-in functionality. Below, we provide a mini dashboard of different quantities and how to generate them. \n", "\n", "If there are certain quantites that you wish to be added that don't exist within FreeGSNKE at the moment, please do submit a feature request to the Github repository." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### General equillibrium quantites\n", "Here, we display a number of used equilibrium quantites. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# plasma quantities\n", "print(rf\"Plasma current: {eq.plasmaCurrent()} [A]\")\n", "print(\"---\")\n", "print(rf\"Poloidal beta (definition 1): {eq.poloidalBeta1()}\")\n", "print(rf\"Poloidal beta (definition 2): {eq.poloidalBeta2()}\")\n", "print(rf\"Poloidal beta (definition 3): {eq.poloidalBeta3()}\")\n", "print(rf\"Poloidal beta (definition 4): {eq.poloidalBeta4()}\")\n", "print(\"---\")\n", "print(rf\"Toroidal beta (definition 1): {eq.toroidalBeta1()}\")\n", "print(rf\"Toroidal beta (definition 2): {eq.toroidalBeta2()}\")\n", "print(rf\"Toroidal beta (definition 3): {eq.toroidalBeta3()}\")\n", "print(rf\"Toroidal beta (definition 4): {eq.toroidalBeta4()}\")\n", "print(\"---\")\n", "print(rf\"Normalised total beta: {eq.normalised_total_Beta()}\")\n", "print(\"---\")\n", "print(fr\"Plasma internal inductance (definition 1): {eq.internalInductance1()}\")\n", "print(fr\"Plasma internal inductance (definition 2): {eq.internalInductance2()}\")\n", "print(fr\"Plasma internal inductance (definition 3): {eq.internalInductance3()}\")\n", "print(fr\"Plasma internal inductance (definition 4): {eq.internalInductance4()}\")\n", "print(fr\"Plasma internal inductance (definition 5): {eq.internalInductance5()}\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# plasma geometry quantities\n", "print(rf\"Minor radius: {eq.minorRadius()} [m]\")\n", "print(\"---\")\n", "print(rf\"Magnetic axis location: {eq.magneticAxis()[0:2]} [m]\")\n", "print(rf\"Geometric axis location: {eq.geometricAxis()} [m]\")\n", "print(rf\"Shafranov shift: {eq.shafranov_shift()} [m]\")\n", "print(rf\"Inner and outer mipdplane (Z=0) radii: {eq.innerOuterSeparatrix(Z=0)} [m]\")\n", "print(rf\"Inner and outer mipdplane (Z=0) radii (alternative): {eq.innerOuterSeparatrix2(Z=0)} [m]\")\n", "\n", "print(\"---\")\n", "print(rf\"LCFS circumference: {eq.separatrix_length()} [m]\")\n", "print(rf\"LCFS area: {eq.separatrix_area()} [m^2]\")\n", "print(rf\"Plasma volume: {eq.plasmaVolume()} [m^3]\")\n", "\n", "print(\"---\")\n", "print(rf\"Aspect ratio: {eq.aspectRatio()}\")\n", "\n", "print(\"---\")\n", "print(rf\"Geometric elongation: {eq.geometricElongation()}\")\n", "print(rf\"Geometric elongation (upper): {eq.geometricElongation_upper()}\")\n", "print(rf\"Geometric elongation (lower): {eq.geometricElongation_lower()}\")\n", "print(rf\"Effective elongation: {eq.effectiveElongation()}\")\n", "\n", "print(\"---\")\n", "print(rf\"Triangularity: {eq.triangularity()}\")\n", "print(rf\"Triangularity (upper): {eq.triangularity_upper()}\")\n", "print(rf\"Triangularity (lower): {eq.triangularity_lower()}\")\n", "\n", "print(\"---\")\n", "s_uo, s_ui, s_lo, s_li = eq.squareness()\n", "print(rf\"Squareness (upper outer): {s_uo}\")\n", "print(rf\"Squareness (upper inner): {s_ui}\")\n", "print(rf\"Squareness (lower outer): {s_lo}\")\n", "print(rf\"Squareness (lower inner): {s_li}\")\n", "\n", "print(\"---\")\n", "L = eq.closest_wall_point()\n", "print(rf\"Point on wall that is closest to the LCFS: {L[0]} [m]\")\n", "print(rf\"Corresponding distance from point to LCFS: {L[1]} [m]\")\n", "\n", "print(\"---\")\n", "print(rf\"Is this a limited plasma?: {eq.flag_limiter}\")\n", "print(rf\"Does the core plasma boundary intersect the wall?: {eq.intersectsWall()}\")\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(fr\"Poloidal flux on magnetic axis: {eq.psi_axis} [Webers/2\\pi]\")\n", "print(fr\"Poloidal flux on plasma boundary: {eq.psi_bndry} [Webers/2\\pi]\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# extract the plasma core boundary (specify the number of points you want)\n", "eq.separatrix(ntheta=20) " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# extract the X points (R, Z, psi)\n", "eq.xpt" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# extract the strike points where the separatric intersects the wall, if any (R, Z)\n", "eq.strikepoints()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# here we plot some of the features above on different axes to show where they are located (see code documentation for more details)\n", "\n", "sep = eq.separatrix(ntheta=360)\n", "strikes = eq.strikepoints()\n", "\n", "fig1, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(15, 8), dpi=80)\n", "plt.subplots_adjust(wspace=0.25) # adjust the horizontal space between subplots\n", "\n", "ax1.grid(True, which='both')\n", "eq.tokamak.plot(axis=ax1,show=False)\n", "ax1.fill(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, facecolor='w', zorder=0)\n", "ax1.contour(eq.R, eq.Z, eq.psi(), levels=[eq.psi_bndry], linestyles='solid', colors='r') #, label=\"Separatrix\")\n", "ax1.scatter(strikes[:,0], strikes[:,1], color='k', marker='x', s=40, zorder=2, label=\"Strikepoints\")\n", "ax1.set_aspect('equal')\n", "ax1.set_xlim(0.1, 2.15)\n", "ax1.set_ylim(-2.25, 2.25)\n", "ax1.set_xlabel(r'Major radius, $R$ [m]')\n", "ax1.set_ylabel(r'Height, $Z$ [m]')\n", "ax1.legend(loc=\"upper right\")\n", "\n", "\n", "\n", "\n", "\n", "ax2.grid(True, which='both')\n", "eq.tokamak.plot(axis=ax2,show=False)\n", "ax2.fill(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, facecolor='w', zorder=0)\n", "ax2.scatter(sep[:,0], sep[:,1], color='k', marker='o', s=2, zorder=2, label=\"LCFS (plasma boundary)\")\n", "ax2.scatter(eq.xpt[:,0], eq.xpt[:,1], color='r', marker='x', s=40, zorder=2, label=\"X-points\")\n", "ax2.scatter(eq.magneticAxis()[0], eq.magneticAxis()[1], color='g', marker='x', s=40, zorder=2, label=\"Magnetic axis\")\n", "ax2.scatter(eq.geometricAxis()[0], eq.geometricAxis()[1], color='b', marker='x', s=40, zorder=2, label=\"Geometric axis\")\n", "ax2.set_aspect('equal')\n", "ax2.set_xlim(0.1, 2.15)\n", "ax2.set_ylim(-2.25, 2.25)\n", "ax2.set_xlabel(r'Major radius, $R$ [m]')\n", "ax2.set_ylabel(r'Height, $Z$ [m]')\n", "ax2.legend(loc=\"upper right\")\n", "\n", "\n", "\n", "\n", "\n", "ax3.grid(True, which='both')\n", "eq.tokamak.plot(axis=ax3,show=False)\n", "ax3.fill(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, facecolor='w', zorder=0)\n", "ax3.scatter(sep[:,0], sep[:,1], color='k', marker='o', s=2, zorder=2, label=\"LCFS (plasma boundary)\")\n", "ax3.scatter(eq._sep_Rmin, eq._sep_ZRmin, color='r', marker='x', s=40, zorder=2, label=\"(Rmin, ZRmin)\")\n", "ax3.scatter(eq._sep_Rmax, eq._sep_ZRmax, color='m', marker='x', s=40, zorder=2, label=\"(Rmax, ZRmax)\")\n", "ax3.scatter(eq._sep_RZmin, eq._sep_Zmin, color='g', marker='x', s=40, zorder=2, label=\"(RZmin, Zmin)\")\n", "ax3.scatter(eq._sep_RZmax, eq._sep_Zmax, color='b', marker='x', s=40, zorder=2, label=\"(ZRmax, Zmax)\")\n", "ax3.set_aspect('equal')\n", "ax3.set_xlim(0.1, 2.15)\n", "ax3.set_ylim(-2.25, 2.25)\n", "ax3.set_xlabel(r'Major radius, $R$ [m]')\n", "ax3.set_ylabel(r'Height, $Z$ [m]')\n", "ax3.legend(loc=\"upper right\")\n", "\n", "plt.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here we visualise the dr sep quantities.\n", "\n", "These two distances are defined as the radial separation (on the inboard and outboard sides) at the vertical position `Z_level` between flux surfaces passing through the lower and upper X-points. The default `Z_level` value is at the geometric axis (i.e. `eq.Zgeometric()`), however, here we have chosen the midplance (at $Z=0$)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "dr_seps = eq.dr_sep(Z_level=0)\n", "\n", "print(f\"dr_sep_in = {dr_seps[0]} [m].\")\n", "print(f\"dr_sep_out = {dr_seps[1]} [m].\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the plot below, dr_sep_in is the difference between the radial position of the blue and the red flux surfaces at $Z=0$. \n", "\n", "Similarly, dr_sep_out is the difference between the radial position of the red and the blue flux surfaces at $Z=0$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import shapely as sh\n", "\n", "# visualising these quantities\n", "psi_bndry = eq.psi_bndry\n", "xpts = eq.xpt\n", "\n", "# define your polygon (example)\n", "polygon = sh.Polygon(np.array([eq.tokamak.wall.R, eq.tokamak.wall.Z]).T)\n", "\n", "# vectorized-ish approach: boolean mask for points inside\n", "mask = np.array([polygon.contains(sh.Point(x, y)) for x, y in xpts[:,0:2]])\n", "\n", "# select points inside as a NumPy array\n", "xpts_inside_wall = xpts[mask,:]\n", "\n", "# get indices of the two cloest to magnetic axis \n", "closest_xpts_idxs = np.argsort(np.linalg.norm(xpts_inside_wall[:,0:2] - eq.opt[0,0:2], axis=1))\n", "closest_xpts_sorted = xpts_inside_wall[closest_xpts_idxs[0:2],:]\n", "\n", "# plot the flux contours for the values of psi_boundary at each X-point\n", "fig1, axes = plt.subplots(1, 3, figsize=(15, 8), dpi=80)\n", "plt.subplots_adjust(wspace=0.25) # adjust the horizontal space between subplots\n", "\n", "axes[0].grid(True, which='both')\n", "eq.tokamak.plot(axis=axes[0],show=False)\n", "axes[0].fill(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, facecolor='w', zorder=0)\n", "axes[0].contour(eq.R, eq.Z, eq.psi(), levels=[closest_xpts_sorted[0, 2]], linestyles='solid', colors='r')\n", "axes[0].contour(eq.R, eq.Z, eq.psi(), levels=[closest_xpts_sorted[1, 2]], linestyles='dashed', colors='b')\n", "axes[0].scatter(closest_xpts_sorted[0,0], closest_xpts_sorted[0,1], color='r', marker='x', s=40, zorder=2, label=\"Separatrix from lower X-point\")\n", "axes[0].scatter(closest_xpts_sorted[1,0], closest_xpts_sorted[1,1], color='b', marker='x', s=40, zorder=2, label=\"Separatrix from upper X-point\")\n", "axes[0].set_aspect('equal')\n", "axes[0].set_xlim(0.1, 2.15)\n", "axes[0].set_ylim(-2.25, 2.25)\n", "axes[0].set_xlabel(r'Major radius, $R$ [m]')\n", "axes[0].set_ylabel(r'Height, $Z$ [m]')\n", "axes[0].legend(loc=\"upper right\")\n", "\n", "\n", "axes[1].grid(True, which='both')\n", "eq.tokamak.plot(axis=axes[1],show=False)\n", "axes[1].fill(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, facecolor='w', zorder=0)\n", "axes[1].contour(eq.R, eq.Z, eq.psi(), levels=[closest_xpts_sorted[0, 2]], linestyles='solid', colors='r')\n", "axes[1].contour(eq.R, eq.Z, eq.psi(), levels=[closest_xpts_sorted[1, 2]], linestyles='dashed', colors='b')\n", "axes[1].scatter(closest_xpts_sorted[0,0], closest_xpts_sorted[0,1], color='r', marker='x', s=40, zorder=2, label=\"Separatrix from lower X-point\")\n", "axes[1].scatter(closest_xpts_sorted[1,0], closest_xpts_sorted[1,1], color='b', marker='x', s=40, zorder=2, label=\"Separatrix from upper X-point\")\n", "axes[1].set_aspect('equal')\n", "axes[1].set_xlim(0.18, 0.38)\n", "axes[1].set_ylim(-0.05, 0.05)\n", "axes[1].set_xlabel(r'Major radius, $R$ [m]')\n", "axes[1].set_ylabel(r'Height, $Z$ [m]')\n", "axes[1].set_title(\"Inboard side: dr_sep_in\")\n", "\n", "\n", "axes[2].grid(True, which='both')\n", "eq.tokamak.plot(axis=axes[2],show=False)\n", "axes[2].fill(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, facecolor='w', zorder=0)\n", "axes[2].contour(eq.R, eq.Z, eq.psi(), levels=[closest_xpts_sorted[0, 2]], linestyles='solid', colors='r')\n", "axes[2].contour(eq.R, eq.Z, eq.psi(), levels=[closest_xpts_sorted[1, 2]], linestyles='dashed', colors='b')\n", "axes[2].scatter(closest_xpts_sorted[0,0], closest_xpts_sorted[0,1], color='r', marker='x', s=40, zorder=2, label=\"Separatrix from lower X-point\")\n", "axes[2].scatter(closest_xpts_sorted[1,0], closest_xpts_sorted[1,1], color='b', marker='x', s=40, zorder=2, label=\"Separatrix from upper X-point\")\n", "axes[2].set_aspect('equal')\n", "axes[2].set_xlim(1.3, 1.5)\n", "axes[2].set_ylim(-0.05, 0.05)\n", "axes[2].set_xlabel(r'Major radius, $R$ [m]')\n", "axes[2].set_ylabel(r'Height, $Z$ [m]')\n", "axes[2].set_title(\"Outboard side: dr_sep_out\")\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Flux quantites\n", "Here, we plot the total poloidal flux $\\psi$ [Webers / $2\\pi$] and its two components $\\psi_p$ (the plasma flux) and $\\psi_c$ (the coil flux): $\\psi = \\psi_p + \\psi_c$. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig1, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(15, 8), dpi=80)\n", "plt.subplots_adjust(wspace=0.5) # adjust the horizontal space between subplots\n", "\n", "ax1.grid(True, which='both')\n", "eq.tokamak.plot(axis=ax1,show=False)\n", "# ax1.plot(tokamak.limiter.R, tokamak.limiter.Z, color='k', linewidth=1.2, linestyle=\"--\")\n", "ax1.plot(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, linestyle=\"-\")\n", "im1 = ax1.contour(eq.R, eq.Z, eq.psi(), levels=50) # total psi\n", "ax1.contour(eq.R, eq.Z, eq.psi(), levels=[eq.psi_bndry], colors='r') # psi boundary contour\n", "ax1.set_aspect('equal')\n", "ax1.set_xlim(0.1, 2.15)\n", "ax1.set_ylim(-2.25, 2.25)\n", "ax1.set_title(r\"Total flux, $\\psi$\")\n", "ax1.set_xlabel(r'Major radius, $R$ [m]')\n", "ax1.set_ylabel(r'Height, $Z$ [m]')\n", "cbar = plt.colorbar(im1, ax=ax1, fraction=0.09)\n", "\n", "\n", "ax2.grid(True, which='both')\n", "eq.tokamak.plot(axis=ax2,show=False)\n", "# ax2.plot(tokamak.limiter.R, tokamak.limiter.Z, color='k', linewidth=1.2, linestyle=\"--\")\n", "ax2.plot(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, linestyle=\"-\")\n", "im2 = ax2.contour(eq.R, eq.Z, eq.plasma_psi, levels=50) # plasma psi\n", "ax2.set_aspect('equal')\n", "ax2.set_xlim(0.1, 2.15)\n", "ax2.set_ylim(-2.25, 2.25)\n", "ax2.set_title(r\"Plasma flux, $\\psi_p$\")\n", "ax2.set_xlabel(r'Major radius, $R$ [m]')\n", "ax2.set_ylabel(r'Height, $Z$ [m]')\n", "cbar = plt.colorbar(im2, ax=ax2, fraction=0.09)\n", "\n", "\n", "ax3.grid(True, which='both')\n", "eq.tokamak.plot(axis=ax3,show=False)\n", "# ax3.plot(tokamak.limiter.R, tokamak.limiter.Z, color='k', linewidth=1.2, linestyle=\"--\")\n", "ax3.plot(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, linestyle=\"-\")\n", "im3 = ax3.contour(eq.R, eq.Z, eq.tokamak_psi, levels=50) # coil psi\n", "ax3.set_aspect('equal')\n", "ax3.set_xlim(0.1, 2.15)\n", "ax3.set_ylim(-2.25, 2.25)\n", "ax3.set_title(r\"Coil flux, $\\psi_c$\")\n", "ax3.set_xlabel(r'Major radius, $R$ [m]')\n", "ax3.set_ylabel(r'Height, $Z$ [m]')\n", "cbar = plt.colorbar(im3, ax=ax3, fraction=0.09)\n", "\n", "\n", "plt.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can also extract the flux produced by individual coils (or pasive structures) if required. We do this my using the pre-calculated Green's functions, the coil currents, multipliers, and polarities. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "currents = eq.tokamak.getCurrents()\n", "\n", "# calculate the active/passive coil fluxes\n", "\n", "psi_coils = dict()\n", "for i, name in enumerate(currents.keys()):\n", " coil = eq.tokamak.coils_dict[name]\n", " scaling = coil[\"multiplier\"]*coil[\"polarity\"]\n", " greens_matrix = 0.0\n", " \n", " if type(eq._pgreen[name]) is dict:\n", " num_coils = len(eq._pgreen[name])\n", " for i, ind in enumerate(eq._pgreen[name]):\n", " greens_matrix += eq._pgreen[name][ind]*scaling[i*(len(scaling)//num_coils)]\n", " else:\n", " num_coils = 1\n", " greens_matrix = eq._pgreen[name]*scaling[0]\n", "\n", " psi_coils[name] = greens_matrix*currents[name]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here, we plot a few of the fluxes from the coils. You can change the `name` parameters to visualise different coil fluxes." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig1, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(15, 8), dpi=80)\n", "plt.subplots_adjust(wspace=0.5) # adjust the horizontal space between subplots\n", "\n", "name = \"Solenoid\"\n", "ax1.grid(True, which='both')\n", "eq.tokamak.plot(axis=ax1,show=False)\n", "# ax1.plot(tokamak.limiter.R, tokamak.limiter.Z, color='k', linewidth=1.2, linestyle=\"--\")\n", "ax1.plot(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, linestyle=\"-\")\n", "im1 = ax1.contour(eq.R, eq.Z, psi_coils[name], levels=50) \n", "ax1.set_aspect('equal')\n", "ax1.set_xlim(0.1, 2.15)\n", "ax1.set_ylim(-2.25, 2.25)\n", "ax1.set_title(rf\"$\\psi_c$ ({name})\")\n", "ax1.set_xlabel(r'Major radius, $R$ [m]')\n", "ax1.set_ylabel(r'Height, $Z$ [m]')\n", "cbar = plt.colorbar(im1, ax=ax1, fraction=0.09)\n", "\n", "name = \"D1\"\n", "ax2.grid(True, which='both')\n", "eq.tokamak.plot(axis=ax2,show=False)\n", "# ax2.plot(tokamak.limiter.R, tokamak.limiter.Z, color='k', linewidth=1.2, linestyle=\"--\")\n", "ax2.plot(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, linestyle=\"-\")\n", "im2 = ax2.contour(eq.R, eq.Z, psi_coils[name], levels=50) \n", "ax2.set_aspect('equal')\n", "ax2.set_xlim(0.1, 2.15)\n", "ax2.set_ylim(-2.25, 2.25)\n", "ax2.set_title(rf\"$\\psi_c$ ({name})\")\n", "ax2.set_xlabel(r'Major radius, $R$ [m]')\n", "ax2.set_ylabel(r'Height, $Z$ [m]')\n", "cbar = plt.colorbar(im2, ax=ax2, fraction=0.09)\n", "\n", "name = \"P6\"\n", "ax3.grid(True, which='both')\n", "eq.tokamak.plot(axis=ax3,show=False)\n", "# ax3.plot(tokamak.limiter.R, tokamak.limiter.Z, color='k', linewidth=1.2, linestyle=\"--\")\n", "ax3.plot(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, linestyle=\"-\")\n", "im3 = ax3.contour(eq.R, eq.Z, psi_coils[name], levels=50) \n", "ax3.set_aspect('equal')\n", "ax3.set_xlim(0.1, 2.15)\n", "ax3.set_ylim(-2.25, 2.25)\n", "ax3.set_title(rf\"$\\psi_c$ ({name})\")\n", "ax3.set_xlabel(r'Major radius, $R$ [m]')\n", "ax3.set_ylabel(r'Height, $Z$ [m]')\n", "cbar = plt.colorbar(im3, ax=ax3, fraction=0.09)\n", "\n", "\n", "plt.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Magnetic fields\n", "\n", "Here, we plot the different magnetic field components over the domain." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig1, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(15, 8), dpi=80)\n", "plt.subplots_adjust(wspace=0.5) # adjust the horizontal space between subplots\n", "\n", "ax1.grid(True, which='both')\n", "eq.tokamak.plot(axis=ax1,show=False)\n", "# ax1.plot(tokamak.limiter.R, tokamak.limiter.Z, color='k', linewidth=1.2, linestyle=\"--\")\n", "ax1.plot(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, linestyle=\"-\")\n", "im1 = ax1.contour(eq.R, eq.Z, eq.Br(eq.R, eq.Z), levels=200) \n", "ax1.set_aspect('equal')\n", "ax1.set_xlim(0.1, 2.15)\n", "ax1.set_ylim(-2.25, 2.25)\n", "ax1.set_title(r\"$B_R$ [T]\")\n", "ax1.set_xlabel(r'Major radius, $R$ [m]')\n", "ax1.set_ylabel(r'Height, $Z$ [m]')\n", "cbar = plt.colorbar(im1, ax=ax1, fraction=0.09)\n", "\n", "\n", "ax2.grid(True, which='both')\n", "eq.tokamak.plot(axis=ax2,show=False)\n", "# ax2.plot(tokamak.limiter.R, tokamak.limiter.Z, color='k', linewidth=1.2, linestyle=\"--\")\n", "ax2.plot(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, linestyle=\"-\")\n", "im2 = ax2.contour(eq.R, eq.Z, eq.Btor(eq.R, eq.Z), levels=100) \n", "ax2.set_aspect('equal')\n", "ax2.set_xlim(0.1, 2.15)\n", "ax2.set_ylim(-2.25, 2.25)\n", "ax2.set_title(r\"$B_{\\phi}$ [T]\")\n", "ax2.set_xlabel(r'Major radius, $R$ [m]')\n", "ax2.set_ylabel(r'Height, $Z$ [m]')\n", "cbar = plt.colorbar(im2, ax=ax2, fraction=0.09)\n", "\n", "\n", "ax3.grid(True, which='both')\n", "eq.tokamak.plot(axis=ax3,show=False)\n", "# ax3.plot(tokamak.limiter.R, tokamak.limiter.Z, color='k', linewidth=1.2, linestyle=\"--\")\n", "ax3.plot(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, linestyle=\"-\")\n", "im3 = ax3.contour(eq.R, eq.Z, eq.Bz(eq.R, eq.Z), levels=100) \n", "ax3.set_aspect('equal')\n", "ax3.set_xlim(0.1, 2.15)\n", "ax3.set_ylim(-2.25, 2.25)\n", "ax3.set_title(r\"$B_Z$ [T]\")\n", "ax3.set_xlabel(r'Major radius, $R$ [m]')\n", "ax3.set_ylabel(r'Height, $Z$ [m]')\n", "cbar = plt.colorbar(im3, ax=ax3, fraction=0.09)\n", "\n", "\n", "plt.tight_layout()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# this is just the contribution from the plasma\n", "fig1, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 8), dpi=80)\n", "plt.subplots_adjust(wspace=0.5) # adjust the horizontal space between subplots\n", "\n", "ax1.grid(True, which='both')\n", "eq.tokamak.plot(axis=ax1,show=False)\n", "# ax1.plot(tokamak.limiter.R, tokamak.limiter.Z, color='k', linewidth=1.2, linestyle=\"--\")\n", "ax1.plot(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, linestyle=\"-\")\n", "im1 = ax1.contour(eq.R, eq.Z, eq.plasmaBr(eq.R, eq.Z), levels=30) \n", "ax1.set_aspect('equal')\n", "ax1.set_xlim(0.1, 2.15)\n", "ax1.set_ylim(-2.25, 2.25)\n", "ax1.set_title(r\"$B_R$ (from plasma) [T]\")\n", "ax1.set_xlabel(r'Major radius, $R$ [m]')\n", "ax1.set_ylabel(r'Height, $Z$ [m]')\n", "cbar = plt.colorbar(im1, ax=ax1, fraction=0.09)\n", "\n", "\n", "ax2.grid(True, which='both')\n", "eq.tokamak.plot(axis=ax2,show=False)\n", "# ax2.plot(tokamak.limiter.R, tokamak.limiter.Z, color='k', linewidth=1.2, linestyle=\"--\")\n", "ax2.plot(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, linestyle=\"-\")\n", "im2 = ax2.contour(eq.R, eq.Z, eq.plasmaBz(eq.R, eq.Z), levels=30) \n", "ax2.set_aspect('equal')\n", "ax2.set_xlim(0.1, 2.15)\n", "ax2.set_ylim(-2.25, 2.25)\n", "ax2.set_title(r\"$B_Z$ (from plasma) [T]\")\n", "ax2.set_xlabel(r'Major radius, $R$ [m]')\n", "ax2.set_ylabel(r'Height, $Z$ [m]')\n", "cbar = plt.colorbar(im2, ax=ax2, fraction=0.09)\n", "\n", "\n", "plt.tight_layout()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# this is just the contribution from the active coils\n", "\n", "Br_actives = 0.0\n", "Bz_actives = 0.0\n", "for name in eq.tokamak.coils_list[0:12]:\n", " Br_actives += eq.tokamak[name].Br(eq.R, eq.Z)\n", " Bz_actives += eq.tokamak[name].Bz(eq.R, eq.Z)\n", "\n", "fig1, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 8), dpi=80)\n", "plt.subplots_adjust(wspace=0.5) # adjust the horizontal space between subplots\n", "\n", "ax1.grid(True, which='both')\n", "eq.tokamak.plot(axis=ax1,show=False)\n", "# ax1.plot(tokamak.limiter.R, tokamak.limiter.Z, color='k', linewidth=1.2, linestyle=\"--\")\n", "ax1.plot(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, linestyle=\"-\")\n", "im1 = ax1.contour(eq.R, eq.Z, Br_actives, levels=200) \n", "ax1.set_aspect('equal')\n", "ax1.set_xlim(0.1, 2.15)\n", "ax1.set_ylim(-2.25, 2.25)\n", "ax1.set_title(r\"$B_R$ (from coils) [T]\")\n", "ax1.set_xlabel(r'Major radius, $R$ [m]')\n", "ax1.set_ylabel(r'Height, $Z$ [m]')\n", "cbar = plt.colorbar(im1, ax=ax1, fraction=0.09)\n", "\n", "\n", "ax2.grid(True, which='both')\n", "eq.tokamak.plot(axis=ax2,show=False)\n", "# ax2.plot(tokamak.limiter.R, tokamak.limiter.Z, color='k', linewidth=1.2, linestyle=\"--\")\n", "ax2.plot(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, linestyle=\"-\")\n", "im2 = ax2.contour(eq.R, eq.Z, Bz_actives, levels=200) \n", "ax2.set_aspect('equal')\n", "ax2.set_xlim(0.1, 2.15)\n", "ax2.set_ylim(-2.25, 2.25)\n", "ax2.set_title(r\"$B_Z$ (from coils) [T]\")\n", "ax2.set_xlabel(r'Major radius, $R$ [m]')\n", "ax2.set_ylabel(r'Height, $Z$ [m]')\n", "cbar = plt.colorbar(im2, ax=ax2, fraction=0.09)\n", "\n", "\n", "plt.tight_layout()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# # this is just the contribution from the passive coils (note this is zero because there are no currents in the passives!)\n", "\n", "# Br_passives = 0.0\n", "# Bz_passives = 0.0\n", "# for name in eq.tokamak.coils_list[12:]:\n", "# Br_passives += eq.tokamak[name].Br(eq.R, eq.Z)\n", "# Bz_passives += eq.tokamak[name].Bz(eq.R, eq.Z)\n", " \n", " \n", "# fig1, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 8), dpi=80)\n", "# plt.subplots_adjust(wspace=0.5) # adjust the horizontal space between subplots\n", "\n", "# ax1.grid(True, which='both')\n", "# eq.tokamak.plot(axis=ax1,show=False)\n", "# # ax1.plot(tokamak.limiter.R, tokamak.limiter.Z, color='k', linewidth=1.2, linestyle=\"--\")\n", "# ax1.plot(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, linestyle=\"-\")\n", "# im1 = ax1.contour(eq.R, eq.Z, Br_passives, levels=200) \n", "# ax1.set_aspect('equal')\n", "# ax1.set_xlim(0.1, 2.15)\n", "# ax1.set_ylim(-2.25, 2.25)\n", "# ax1.set_title(r\"$B_R$ (from coils) [T]\")\n", "# ax1.set_xlabel(r'Major radius, $R$ [m]')\n", "# ax1.set_ylabel(r'Height, $Z$ [m]')\n", "# cbar = plt.colorbar(im1, ax=ax1, fraction=0.09)\n", "\n", "\n", "# ax2.grid(True, which='both')\n", "# eq.tokamak.plot(axis=ax2,show=False)\n", "# # ax2.plot(tokamak.limiter.R, tokamak.limiter.Z, color='k', linewidth=1.2, linestyle=\"--\")\n", "# ax2.plot(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, linestyle=\"-\")\n", "# im2 = ax2.contour(eq.R, eq.Z, Bz_passives, levels=200) \n", "# ax2.set_aspect('equal')\n", "# ax2.set_xlim(0.1, 2.15)\n", "# ax2.set_ylim(-2.25, 2.25)\n", "# ax2.set_title(r\"$B_Z$ (from coils) [T]\")\n", "# ax2.set_xlabel(r'Major radius, $R$ [m]')\n", "# ax2.set_ylabel(r'Height, $Z$ [m]')\n", "# cbar = plt.colorbar(im2, ax=ax2, fraction=0.09)\n", "\n", "\n", "# plt.tight_layout()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# total poloidal magnetic field\n", "\n", "fig1, ax1 = plt.subplots(1, 1, figsize=(4, 8), dpi=100)\n", "\n", "ax1.grid(True, which='both')\n", "eq.tokamak.plot(axis=ax1,show=False)\n", "# ax1.plot(tokamak.limiter.R, tokamak.limiter.Z, color='k', linewidth=1.2, linestyle=\"--\")\n", "ax1.plot(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, linestyle=\"-\")\n", "im1 = ax1.contour(eq.R, eq.Z, eq.Bpol(eq.R, eq.Z), levels=200) \n", "ax1.set_aspect('equal')\n", "ax1.set_xlim(0.1, 2.15)\n", "ax1.set_ylim(-2.25, 2.25)\n", "ax1.set_title(r\"$B_{pol}$ [T]\")\n", "ax1.set_xlabel(r'Major radius, $R$ [m]')\n", "ax1.set_ylabel(r'Height, $Z$ [m]')\n", "cbar = plt.colorbar(im1, ax=ax1, fraction=0.09)\n", "\n", "plt.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Coil currents\n", "\n", "Here, we visualise the size of the currents in the poloidal field coils." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# listing the currents may not be that informative \n", "eq.tokamak.getCurrents()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You can play around with the colour schemes below to visualise how the current is distributed around the machine. Note here that there is no current in the passive structures but they can be included if present. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from matplotlib.patches import Rectangle, Polygon\n", "from matplotlib.colors import Normalize, TwoSlopeNorm\n", "import matplotlib.cm as cm\n", "\n", "\n", "\n", "# create colormap based on magnitude of currents\n", "currents_array = []\n", "for key in list(eq.tokamak.coils_dict.keys()):\n", " currents_array.append(eq.tokamak[key].current)\n", "\n", "max_curr = np.max(np.abs(currents_array))\n", "# norm = Normalize(vmin=-max_curr, vmax=max_curr) # alternative colorbar\n", "norm = TwoSlopeNorm(vmin=np.min(currents_array), vcenter=0, vmax=np.max(currents_array)) \n", "cmap = cm.bwr\n", "\n", "\n", "# plot\n", "fig1, ax1 = plt.subplots(1, 1, figsize=(4, 8), dpi=120)\n", "plt.tight_layout()\n", "ax1.plot(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, linestyle=\"-\")\n", "ax1.grid(alpha=0.5)\n", "ax1.set_aspect('equal')\n", "ax1.set_xlim(0.1, 2.15)\n", "ax1.set_ylim(-2.25, 2.25)\n", "ax1.set_xlabel(r'Major radius, $R$ [m]')\n", "ax1.set_ylabel(r'Height, $Z$ [m]')\n", "\n", "\n", "for name in list(eq.tokamak.coils_dict.keys()):\n", " coil = eq.tokamak.coils_dict[name]\n", " current = eq.tokamak[name].current\n", " color = cmap(norm(current)) # map the current to a color\n", " \n", " # plot active coils (and currents)\n", " if coil[\"active\"]:\n", "\n", " for i in range(0, len(coil[\"coords\"][0,:])):\n", " patch = Rectangle(\n", " (coil[\"coords\"][0,i] - coil[\"dR\"] / 2, coil[\"coords\"][1,i] - coil[\"dZ\"] / 2),\n", " width=coil[\"dR\"],\n", " height=coil[\"dZ\"],\n", " facecolor=color,\n", " edgecolor='k',\n", " linewidth=0.2,\n", " )\n", " ax1.add_patch(patch)\n", " \n", " # plot passive structures (currents are zero here) \n", " else:\n", " patch = Polygon(\n", " coil[\"vertices\"].T, \n", " facecolor=color, \n", " edgecolor=\"k\", \n", " linewidth=0.2\n", " )\n", " ax1.add_patch(patch)\n", "\n", "\n", "# add a colorbar\n", "sm = cm.ScalarMappable(cmap=cmap, norm=norm)\n", "sm.set_array([]) \n", "cbar = fig1.colorbar(sm, ax=ax1, orientation='vertical', fraction=0.09)\n", "cbar.set_label('Current [A]')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Forces on the coils\n", "\n", "We can also extract the radial and vertical forces on the coils." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "eq.printForces()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### 1D plasma current density profiles (and others)\n", "Here, we visualise the $p'$ and $FF'$ profiles used in our equilbirium solve. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# plot the p' and FF' profiles\n", "\n", "psi_n = eq.psiN_1D(N=65)\n", "\n", "fig1, (ax1, ax2) = plt.subplots(1, 2, figsize=(15,6), dpi=80)\n", "ax1.grid(zorder=0, alpha=0.75)\n", "ax1.plot(psi_n, profiles.pprime(psi_n), color='k', linewidth=1, marker='x', markersize=2, zorder=10)\n", "ax1.set_xlabel(r'$\\hat{\\psi}$')\n", "ax1.set_ylabel(r\"$p'(\\hat{\\psi})$\")\n", "ax1.ticklabel_format(axis='y', scilimits=(0,0))\n", "\n", "ax2.grid(zorder=0, alpha=0.75)\n", "ax2.plot(psi_n, profiles.ffprime(psi_n), color='k', linewidth=1, marker='x', markersize=2, zorder=10)\n", "ax2.set_xlabel(r'$\\hat{\\psi}$')\n", "ax2.set_ylabel(r\"$FF'(\\hat{\\psi})$\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# plot the p and F profiles\n", "\n", "psi_n = eq.psiN_1D(N=65)\n", "\n", "fig1, (ax1, ax2) = plt.subplots(1, 2, figsize=(15,6), dpi=80)\n", "ax1.grid(zorder=0, alpha=0.75)\n", "ax1.plot(psi_n, profiles.pressure(psi_n), color='k', linewidth=1, marker='x', markersize=2, zorder=10)\n", "ax1.set_xlabel(r'$\\hat{\\psi}$')\n", "ax1.set_ylabel(r\"$p(\\hat{\\psi})$\")\n", "ax1.ticklabel_format(axis='y', scilimits=(0,0))\n", "\n", "ax2.grid(zorder=0, alpha=0.75)\n", "ax2.plot(psi_n, profiles.fpol(psi_n), color='k', linewidth=1, marker='x', markersize=2, zorder=10)\n", "ax2.set_xlabel(r'$\\hat{\\psi}$')\n", "ax2.set_ylabel(r\"$F(\\hat{\\psi})$\")\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# plot the p and F profiles\n", "\n", "psi_n = eq.psiN_1D(N=65)\n", "\n", "fig1, (ax1, ax2) = plt.subplots(1, 2, figsize=(15,6), dpi=80)\n", "ax1.grid(zorder=0, alpha=0.75)\n", "ax1.plot(psi_n, profiles.pressure(psi_n), color='k', linewidth=1, marker='x', markersize=2, zorder=10)\n", "ax1.set_xlabel(r'$\\hat{\\psi}$')\n", "ax1.set_ylabel(r\"$p(\\hat{\\psi})$\")\n", "ax1.ticklabel_format(axis='y', scilimits=(0,0))\n", "\n", "ax2.grid(zorder=0, alpha=0.75)\n", "ax2.plot(psi_n, profiles.fpol(psi_n), color='k', linewidth=1, marker='x', markersize=2, zorder=10)\n", "ax2.set_xlabel(r'$\\hat{\\psi}$')\n", "ax2.set_ylabel(r\"$F(\\hat{\\psi})$\")\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# plot q profile\n", "\n", "psi_n = np.linspace(0.01,0.99,65) # values of q at 0 and 1 can be problematic\n", "\n", "fig1, ax1 = plt.subplots(1, 1, figsize=(6,6), dpi=80)\n", "ax1.grid(zorder=0, alpha=0.75)\n", "ax1.plot(psi_n, eq.q(psi_n), color='k', linewidth=1, marker='x', markersize=2, zorder=10)\n", "ax1.set_xlabel(r'$\\hat{\\psi}$')\n", "ax1.set_ylabel(r\"$q(\\hat{\\psi})$\")\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Toroidal current density\n", "\n", "Here, we visualise the toroidal current density inside the core of the plasma. \n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig1, ax1 = plt.subplots(1, 1, figsize=(4, 8), dpi=100)\n", "\n", "ax1.grid(True, which='both')\n", "eq.tokamak.plot(axis=ax1,show=False)\n", "ax1.plot(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, linestyle=\"-\")\n", "im1 = ax1.contourf(eq.R, eq.Z, eq._profiles.jtor, levels=np.linspace(0.01, np.max( eq._profiles.jtor,), 20)) \n", "ax1.contour(eq.R, eq.Z, eq.psi(), levels=[eq.psi_bndry], colors='r')\n", "ax1.set_aspect('equal')\n", "ax1.set_xlim(0.1, 2.15)\n", "ax1.set_ylim(-2.25, 2.25)\n", "ax1.set_title(r\"Toroidal current density, $J_p$ [A]\")\n", "ax1.set_xlabel(r'Major radius, $R$ [m]')\n", "ax1.set_ylabel(r'Height, $Z$ [m]')\n", "cbar = plt.colorbar(im1, ax=ax1, fraction=0.09)\n", "\n", "plt.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Masking arrays\n", "There are a number of masking arrays that are built during the equilibrium solve that may be useful for plotting or other calculations you may carry out. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig1, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 8), dpi=80)\n", "plt.subplots_adjust(wspace=0.25) # adjust the horizontal space between subplots\n", "\n", "ax1.grid(True, which='both')\n", "# eq.tokamak.plot(axis=ax1,show=False)\n", "# ax1.fill(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, facecolor=None, zorder=0)\n", "# ax1.contour(eq.R, eq.Z, eq.psi(), levels=[eq.psi_bndry], linestyles='solid', colors='r') #, label=\"Separatrix\")\n", "ax1.pcolormesh(eq.R, eq.Z, eq.mask_inside_limiter, cmap=\"Greys\")\n", "ax1.set_aspect('equal')\n", "ax1.set_xlim(0.1, 2.15)\n", "ax1.set_ylim(-2.25, 2.25)\n", "ax1.set_xlabel(r'Major radius, $R$ [m]')\n", "ax1.set_ylabel(r'Height, $Z$ [m]')\n", "ax1.set_title(\"Mask of the region allowed to the plasma\")\n", "\n", "\n", "ax2.grid(True, which='both')\n", "# eq.tokamak.plot(axis=ax2,show=False)\n", "# ax2.fill(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, facecolor=None, zorder=0)\n", "# ax2.contour(eq.R, eq.Z, eq.psi(), levels=[eq.psi_bndry], linestyles='solid', colors='r') #, label=\"Separatrix\")\n", "ax2.pcolormesh(eq.R, eq.Z, eq._profiles.diverted_core_mask, cmap=\"Greys\")\n", "ax2.set_aspect('equal')\n", "ax2.set_xlim(0.1, 2.15)\n", "ax2.set_ylim(-2.25, 2.25)\n", "ax2.set_xlabel(r'Major radius, $R$ [m]')\n", "ax2.set_ylabel(r'Height, $Z$ [m]')\n", "ax2.set_title(\"Mask within the plasma core boundary\")\n", "\n", "plt.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Flux averaged quantities\n", "\n", "The `equilibrium` class provides a method to calculate the \"flux averaged\" value of a user-defined 2D scalar field $f(R,Z)$, using [line integrals](https://tutorial.math.lamar.edu/classes/calciii/LineIntegralsPtI.aspx), on a given (normalised) flux surface of $\\psi_n$ (within the last closed flux surface). The flux average $\\langle f \\rangle$ is given by\n", "\n", "$$\n", "\\langle f \\rangle (\\psi_n) = \\frac{ \\int_{C(\\psi_n)} \\frac{f(R,Z)}{B_{\\text{pol}}(R,Z)}\\, ds}{ \\int_{C(\\psi_n)} \\frac{1}{B_{\\text{pol}}(R,Z)} \\, ds },\n", "$$\n", "\n", "where:\n", "- $f(R,Z)$ = 2D scalar field function (e.g. the plasma current density function $J_p(R,Z)$).\n", "- $\\psi_n$ = value of normalised flux at which to evaluate line integrals.\n", "- $C(\\psi_n)$ = curve of $(R, Z)$ points satisfying $\\psi_n = \\text{const}$ (i.e. a flux contour).\n", "- $B_{\\text{pol}}(R,Z)$ = 2D scalar poloidal magnetic field function.\n", "- $ds$ = the arc length element (where $ds = \\sqrt{(dR/dl)^2 + (dZ/dl)^2} dl$ and $l \\in [0,L]$ is a parameterised length going from the beginning to the end of the contour).\n", "\n", "The definition of the flux average was taken from [Song et al. (2024)](https://www.mdpi.com/2571-6182/7/4/45)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here we will calculate the flux averaged values of the plasma current density by first defining a function that returns the current density at arbitrary $(R,Z)$ locations. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from scipy.interpolate import RectBivariateSpline\n", "\n", "def f(R,Z):\n", " jtor = RectBivariateSpline(eq.R_1D, eq.Z_1D, eq._profiles.jtor)\n", " return jtor(R, Z, grid=False)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next we can call the method in the equilibrium object and plot the results. Given the notation above, we note that $\\psi_n$ and $\\hat{\\psi}$ are equivalent. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# call the method\n", "flux_averaged_jtor, psi_n = eq.flux_averaged_function(\n", " f=f,\n", " psi_n=np.linspace(0.0,1.0,101)\n", " )" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# plot\n", "fig1, ax1 = plt.subplots(1, 1, figsize=(6,6), dpi=80)\n", "ax1.grid(zorder=0, alpha=0.75)\n", "ax1.plot(psi_n, flux_averaged_jtor, color='k', linewidth=1, marker='x', markersize=2, zorder=10)\n", "ax1.set_xlabel(r'$\\hat{\\psi}$')\n", "ax1.set_ylabel(r\"$\\langle J_p \\rangle (\\hat{\\psi})$\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# # for example you could flux average other quantities of interest\n", "\n", "# # 1/R\n", "# def f(R,Z):\n", "# g = RectBivariateSpline(eq.R_1D, eq.Z_1D, 1/eq.R)\n", "# return g(R, Z, grid=False)" ] } ], "metadata": { "kernelspec": { "display_name": "freegsnke", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10.16" } }, "nbformat": 4, "nbformat_minor": 4 }